
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c    routine to invert a matrix by Gaussian elimination
c     Ainv=inverse(A)    A(n,n)   Ainv(n,n)
c     Note: sizes maxed at 64 x 64 - re dimension for larger matrices;
c           D matrix is extended.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine invert( A, Ainv, n )
       implicit double precision (a-h,o-z)
       double precision A(n,n),Ainv(n,n)
       double precision D(n,2*n)


c  initialize the reduction matrix
       n2 = 2*n
       do  i = 1,n
        do  j = 1,n
         D(i,j) = A(i,j)
         d(i,n+j) = 0.
        enddo
        D(i,n+i) = 1.
       enddo

c  do the reduction 
       do 3 i = 1,n
        alpha = D(i,i)
        if(alpha .eq. 0.) go to 300
        do 4 j = 1,n2
         D(i,j) = D(i,j)/alpha
 4      continue
        do 5 k = 1,n
         if((k-i).eq.0) go to 5
         beta = D(k,i)
         do 6 j = 1,n2
          D(k,j) = D(k,j) - beta*D(i,j)
 6       continue
 5      continue
 3     continue

c  copy result into output matrix
       do 7 i = 1,n
        do 8 j = 1,n
         Ainv(i,j) = D(i,j+n)
 8      continue
 7     continue


       return
 300   print *,'*** ERROR: Singular matrix ***'
       Ainv=0.0d0 !added by Xiaodong
       return
       end
